/*Filipe Augusto de Souza Fragoso  Nusp - 7277454*/

/*ALGORITMO PARALELO DE GAUSS LENGENDRE PARA APROXIMACAO DO VALOR DO PI*/

/* Execucao do programa:  gcc G-P.c -o G-P.exe -lgmp -pthread 
   Para saber quantas casas de precisao foi obtido executar: ./casa.exe */

#include <stdio.h>
#include <stdlib.h>
#include <pthread.h>
#include <time.h>
#include <gmp.h>
#define DIGITOS 10000000   //numero de casas do pi que sera salvo
#define IT 23   //numero de iteracoes do algoritmo 

/* declaracao de variaveis GMP
Nro IT + 1 de casas devido a (inicializacao + nro de iteracoes )*/
mpf_t a[IT+1],b[IT+1],t[IT+1],p[IT+1],pi[IT+1];
/* variaveis auxiliares */
mpf_t aux1,aux2,auxt1,auxt2,auxt4_1,auxt4_2;

// funcao da thread corresponde a t[i+1]
void * exec_thread1(void *id)
{
	
	int *id_ptr, i;
   	id_ptr = (int *) id;
   	i = *id_ptr;
	
	mpf_sub(aux1,a[i],a[i+1]);      //t[i+1] = t[i] - ( p[i] * ( pow( a[i] - a[i+1], 2 ) ) );
      	mpf_pow_ui (aux1,aux1,2);	
      	mpf_mul(aux2,p[i],aux1);        
	mpf_sub(t[i+1],t[i],aux2); 

	
	return NULL;
}

// funcao da thread corresponde a b[i+1]
void * exec_thread2(void *id)
{
	int *id_ptr, i;
   	id_ptr = (int *) id;
   	i = *id_ptr;

	mpf_mul(auxt2,a[i],b[i]);        //b[i+1] = sqrt( a[i] * b[i]);
      	mpf_sqrt(b[i+1], auxt2);   	
	
	return NULL;
}

// funcao da thread corresponde a p[i+1]
void * exec_thread3(void *id)
{
	int *id_ptr, i;
   	id_ptr = (int *) id;
   	i = *id_ptr;

   	mpf_mul_ui(p[i+1],p[i],2);      //p[i+1] = 2 * p[i];     	
	
	return NULL;
}

// funcao da thread corresponde a pi[i+1]
void * exec_thread4(void *id)
{
	int *id_ptr, i;
   	id_ptr = (int *) id;
   	i = *id_ptr;

      	mpf_add(auxt4_1,a[i],b[i]);   
      	mpf_pow_ui (auxt4_1,auxt4_1,2);       // pi[i] =  pow(a[i]+ b[i],2)/(4 * t[i]);
      	mpf_mul_ui(auxt4_2,t [i],4); 
      	mpf_div(pi[i+1],auxt4_1,auxt4_2);    
	
	return NULL;
}



int main(int argc, char *argv[])
{	
	/* variavel usada para poder ser passado a iteracao atual às funcoes das threads */
	int *taskids[IT+1];
		
	/* declaracao de threads*/
	pthread_t thread[5];
    	
	/* variavel usada para calcular o tempo inicial do programa*/	
	time_t ini = time(0);		
    	
	/* variavel usada para iteracoes*/
	int i = 0;

	/* setando a precisao*/
	//mpf_set_default_prec((DIGITOS*4)+16);
	mpf_set_default_prec(33300000);
        /* inicializacao das variaveis auxiliares do tipo GMP*/	
	mpf_init(aux1);
	mpf_init(aux2);	
	mpf_init(auxt1);
	mpf_init(auxt2);
	mpf_init(auxt4_1);
	mpf_init(auxt4_2);

	/* declaracao dos valores inicias do algoritmo */

	mpf_init_set_ui(a[0], 1);      //valor inicial de a = 1

	mpf_init_set_ui(b[0], 1);      //valor inicial de b = 1/ sqrt(2)
	mpf_sqrt_ui(aux1, 2);   	
	mpf_div(b[0],b[0],aux1);	

	mpf_init_set_d(t[0], 0.25);   //valor inicial de t = 0.25
	
	mpf_init_set_ui(p[0], 1);      //valor inicial de p = 1

	mpf_init_set_ui(pi[0], 0);     //valor inicial de pi = 0


    /*ITERACOES*/
    while(i<IT)
    {
        /* inicializacao das variaveis do tipo GMP na posicao i+1 dos vetores*/
	mpf_init(a[i+1]);
	mpf_init(b[i+1]);	
	mpf_init(t[i+1]);
	mpf_init(p[i+1]);
	mpf_init(pi[i+1]);
	
	// primeiramente calcula-se o valor de a[i+1] pois t[i+1] necessita desse valor para seu calculo. 
	mpf_add(auxt1,a[i],b[i]);
      	mpf_div_ui(auxt1, auxt1,2);       //a[i+1] = ( a[i] + b[i] ) / 2;
      	mpf_set(a[i+1], auxt1);

	taskids[i] = (int *) malloc(sizeof(int));
  	*taskids[i] = i;

	/*Como b[i+1],p[i+1],pi[i+1] necessitam somente de valores da iteracao anterior, e o valor de a[i+1] ja foi calculado,
	  podemos executar b[i+1],p[i+1],pi[i+1] e t[i+1] em paralelo atraves de threads */

	/* cria as threads */	
	pthread_create(&thread[1], NULL, exec_thread1,(void *) taskids[i]);  //thread corresponde a t[i+1]      
	pthread_create(&thread[2], NULL, exec_thread2,(void *) taskids[i]);  //thread corresponde a b[i+1]
	pthread_create(&thread[3], NULL, exec_thread3,(void *) taskids[i]);  //thread corresponde a p[i+1]
	pthread_create(&thread[4], NULL, exec_thread4,(void *) taskids[i]);  //thread corresponde a pi[i+1]
 
	/* espera as threads terminarem */		
	pthread_join(thread[1], NULL);
	pthread_join(thread[2], NULL);
	pthread_join(thread[3], NULL);
	pthread_join(thread[4], NULL);	
	

	i++;     
    }
	/* calcula o tempo de execucao do algoritmo */	
	int tempo_d = time(0)-ini;  	
	/* imprime o tempo de execucao do algoritmo */
	printf("tempo de execucao:%dsegs \n",tempo_d);
	/* tempo de execucao eh calculado antes de escrever o valor de pi no arquivo para saber-se 
	somente o tempo de execucao do algoritmo de Gauss Lengendre*/

	/* ponteiro para o arquivo*/
	FILE *fp;	

	fp = fopen ("exit_p.txt", "w");
	
	/* verifica se arquivo foi aberto com sucesso */
	if(fp == NULL)
	{
		printf("\n Erro de arquivo \n");
	}

	else 
	{
	/* escreve o valor de pi com 10milhoes de casas de precisao no arquivo apontado por fp */
	mpf_out_str(fp,10,33300000, pi[i]);	
	fclose(fp);
	}

			   
  return 0;
}

